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Abstract 

A method for symbolically computing conservation laws of nonlinear partial differential 
equations (PDEs) in multiple space dimensions is presented in the language of variational 
calculus and linear algebra. The steps of the method are illustrated using the Zakharov- 
Kuznetsov and Kadomtsev-Petviashvili equations as examples. 

1-1-1 m 

The method is algorithmic and has been implemented in Mathematica. The software 
package, ConservationLawsMD.m, can be used to symbolically compute and test con- 
servation laws for polynomial PDEs that can be written as nonlinear evolution equations. 

The code ConservationLawsMD.m has been applied to multi-dimensional versions of 
the Sawada-Kotera, Camassa-Holm, Gardner, and Khokhlov-Zabolotskaya equations. 
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Many nonlinear partial differential equations (PDEs) in the applied sciences and engi- 
neering are continuity equations which express conservation of mass, momentum, energy, or 
electric charge. Such equations occur in, e.g., fluid mechanics, particle and quantum physics, 
^ plasma physics, elasticity, gas dynamics, electromagnetism, magneto-hydro-dynamics, non- 

linear optics, etc. Certain nonlinear PDEs admit infinitely many conservation laws. Al- 
though most lack a physical interpretation, these conservation laws play an important role 
in establishing the complete integrability of the PDE. Completely integrable PDEs are nonlin- 
ear PDEs that can be linearized by some transformation (e.g., the Cole-Hopf transformation 
linearizes the B urgers equation) or expl i citly s olved with the Inverse Scattering Transform 
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(1ST). See, e.g.. lAblowitz and Clarksonl (J1991I ). 



The search for conservation laws of the Korteweg-de Vries (KdV) equation began around 
1964 and the knowledge of conservation laws was paramount for the development of soliton 
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theory. As iNewelll ( 119831 ) narrates, the study of conservation laws led to the discovery of 
the Miura transformation (wh ich connect s solutions of the KdV and modified KdV (mKdV) 
equations) and the Lax pair ( Lax! . Il968l ). i.e., a system of linear equations which are only 
compatible i f the original nonlinear PDE holds. In turn, the Lax pai r is the starting point 
for the 1ST (jAblowitz and Clarksonl . Il99ll ; lAblowitz and Segurl . Il98lh which has been used 
to construct soliton solutions, i.e., stable solutions that interact elastically upon collision. 

Conversely, the existence of many (independent) conserved densities is a predictor for 
complete integrability. The knowledge of conservation laws also aids the study of qualita- 
tive properties of PDEs, in p articular, bi-Hamiltonian structures and recursion operators 
( Baldwin and Heremanl . |2010| ). Furthermore, if constitutive properties have been added to 
close "a model," one should verify that conserve d quantities have r emained intact. Another 
application involves numerical solvers for PDEs (ISanz-Sernal . Il982l ). where one checks if the 
first few (discretized) conserved densities are preserved after each time step. 

There are s everal methods for computin g cons e rvation laws as d i scusse d by e.g., 

Bluman et all (l2010h . IJiereman et all d2005h . lNazl (l2008t ). lNaz et all (J2008h . and iRosenhaus 



(120021 ). One could apply Noether's theorem, which states that a (variational) symmetry of 
the PDE corresponds to a conservation law. Using Noether's method, the Differential- 
Geometry package in Maple contains t o ols fo r conservation laws developed by lAnderson 



Wolf 



(2004b) and lAnderson and Cheb-Terrabl (120091 ) . Circumventing Noether's theorem, 
(120021 ) has developed four programs in REDUCE which solve an over-determined system 
of different i al equation s to get conservation laws. Based on the integrating factor method, 



Cheviakovi (120071 . I2OIO1 ) has written a Maple program that computes a set of integrating fac- 



tors (multipliers) on the PDE. To find conservation laws, here again, one has to solve a sys- 
tem o f differential equations. The Maple package PDEtools by ICheb-Terrab and von Bulow 
(120041 ) has the commands ConservedCurrents and ConservedCurrentTest for com- 
puting and testing conservation laws using the integrating fact or method. Last, conservation 



Zakharov and Shabatl (119721 ) 



law s can be obtained from t he La x operators, as shown by, e.g. 
and lDrinfel'd and Sokolov! (119851 ). 

By contrast, the method discussed in this paper uses tools from calculus, the calculus of 
variations, linear algebra, and differential geometry. Briefly, our method works as follows. 
A candidate (local) density is assumed to be a linear combination with undetermined co- 
efficients of monomials that are invariant under the scaling symmetry of the PDE. Next, 
the time derivative of the candidate density is computed and evaluated on the PDE. Sub- 
sequently, the variational derivative is applied to get a linear system for the undetermined 
coefficients. The solution of that system is substituted into the candidate density. Once 
the density is known, the flux is obtained by applying a homotopy operator to invert a 
divergence. Our method can be implemented in any major comput e r alge bra system (CAS). 
The package CqnservationLaw sMD.m bvlPoole and Heremanl (120091 ) is a Mathe matica 
i mplem ent at ion based on work by lHereman et all ( 20051 ) . with new features added by IPoold 
(l2009h . 

This paper is organized as follows. To set the stage, Section [2] shows conservation laws 
for the Zakharov-Kuznetsov (ZK) and Kadomtsev-Petviashvili (KP) equations. Section [3] 
covers the tools that will be used in the algorithm. In Section HJ the algorithm is presented 
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and illustrated for the ZK and KP equations. Section [6] discusses conservation laws of 
PDEs in multiple space dimensions, including the Khokhlov-Zabolotskaya (KZ) equation and 
multi-dimensional versions of the Sawada-Kotera (SK), Camassa-Holm (CH) and Gardner 
equations. Conservation laws for the mult i- dimensional SK, CH, and Gardner equations 
were not found in a literature survey and are presented here for the first time. A general 
conservation law for the KP equation is given in Section Using the (2+l)-dimensional 
Gardner equation as an example, Section [7] shows how to use ConservationLawsMD.m. 
Finally, some conclusions are drawn in Section [HJ 



2. Examples of Conservation Laws 

This paper deals with systems of polynomial PDEs of order M, 

A(u( M )(x)) = 0, 



in n dimensions where x = (x 1 , x 2 , . . . , x n ) is the independent variable. u^ M ^(x) denotes the 
dependent variable u = (u 1 , . . . , v? , . . . , u N ) and its partial derivatives (up to order M) with 
respect to x. We do not cover systems of PDEs with variable coefficients. 
A conservation law for ([!]) is a scalar PDE in the form 



Div P = on A = 0, 



(2) 



where P = P(x, u^ p ) (x)) of some order P. The definition follows lOlyeij ( 19931 ) and lBluman et al 



(120 lOl ). and is commonly used in literature on symmetries of PDEs. In physics, P is called 
a conserved current. More precisely, a conse rvation law can be viewed as an equivalence 
class of conserved currents (jVinogradovl . Il989l ). Our algorithm computes one member from 
each equivalence class; usually a representative that is of lowest complexity and free of curl 
terms. 

Since we work on PDEs from the physical sciences, the algorithm and code are restricted 
to 1 D, 2 D, and 3 D in space, but can be extended to n dimensions. Indeed, many of our 
applications model dynamical problems, where x = (x, y, t) for PDEs in 2 D or x = (x, y, z, t) 
for PDEs in 3 D in space. In either case, the additional variable, t, denotes time. 

Throughout the paper, we will use an alternative definition for ([2]), 



£> 4 p + DivJ = onA = 0, 



(3) 



where p = p(x, u^^(x)) is the conserve d density of s ome order K, and J = J (x, u^Hx) ) 
is the associated flux of some order L ( iMiura et all Il968l ; lAblowitz and Clarksonl . Il99l[ ). 
Comparing §Z§ with <^j, it should be clear that P = (p, J) with P = max{K, L}. 

For simplicity, in the examples we will denote the dependent variables -u 1 ,?/ 2 ,^, etc., 
by u,v,w, etc. Partial derivatives are denoted by subscripts, e.g., f xk \ d \ 2 dzH ^ s wr itten 
as Uk lX k 2 yk 3 z, where the hi are non-negative integers. In ([3]), Div J is the total divergence 
operator, where Div J = V X J X + V y jv if J = (J x , jv) and Div J = V X J X + V y .P + V Z J Z if 
J = (J x , J y , J z ). Logically, D t , T> x , V y , and T> z are total derivative operators. For example, 



the total derivative operator V x (in ID) acting on / = f(x,t,u^ M \x : t)) of order M is 
defined as 

j=l fc=0 OU kx 

where M\ is the order of / in component u^ and M = maxjM-j 1 , . . . , Mf}. The partial 
derivative J^ acts on any x that appears explicitly in /, but not on vP or any partial 
derivatives of vP. Total derivative operators in multiple dimensions are defined analogously 
(see Section |3J). 

The algorithm described in Section H] allows one to compute local conservation laws for 
systems of nonlinear PDEs that can be written as evolution equations. For example, if 
x = (x, y, z, t), an evolution equation in variable t has the form 

u t = G(it ,u x ,u y ,u z ,u 2x ,U2y,u 2z ,u xy , . . . ,u M N xM N yM N z j, (5) 

where G is assumed to be smooth and M\ , M 3 2 , and M| are the orders of component u^ with 
respect to x, y, and z, respectively, and M is the maximum total order of all terms in the 
differential function. Few mult i- dimensional systems of PDEs are of the form (J5J). However, 
it is often possible to obtain a systems of evolution equations by recasting a single higher- 
order equation into a system of first-order equations, sometimes in conjunction with a simple 
transformation. If necessary, our program internally interchanges independent variables to 
obtain (J5J), where time is the evolution variable. However, that swap of variables is not 
used in this paper. For a clearer description of the algorithm, we allow systems of evolution 
equations where any component of x can play the role of evolution variable. 

We now introduce two well-documented PDEs together with some of their conservation 
laws. These PDEs will be used in Section H] to illustrate the steps of the algorithm. 

Example 1. The Zakharov-Kuznetsov (ZK) equation is an evolution equation that models 
three-dimensional ion-sound solitons in a low pressure uniform magnetized plasma 



(IZakharov and Kuznetsovl . 11974 ). After re-scaling, it takes the form 

u t + auu x + (3(Au) x = 0, (6) 

where w(x) = u(x,y, z,t),a and (3 are real parameters, and A = -j^ + jpz + Jp- is the 
Laplacian in 3D. The conservation laws for the (2+l)-dimensional ZK equation, 

Ut + auu x + (3(u 2x + u 2y ) x = 0, (7) 



whe re u(x) = u(x, y, t ) , were studied by, e.g.. lZakharov and Kuznetsovl (11974 ulnf eld (1985 



andlShivamoggi et a/.l ( 1l993l ). After correcting some of the results reported in 



Shivamoggi et all 



( 119931 ). the polynomial conservation laws of (J7|) are 

V t (u) + V x (\au 2 + (3u 2x ) + V y ((3u xy ) = 0, 

which corresponds to the ZK equation itself, and 
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V t (V) +V X (jcm 3 - (3(u 2 x - u 2 y ) + 2(3u(u 2x + u 2y )^ +V y (- 2(3u x u y ^ = 0, (9) 

V t (u 3 - 3f (t£ + ul)) +V x (?>u 2 {\au 2 + (3u 2x ) - 6(3u(u 2 x + u 2 y ) + 3^(u 2 2x - u\ y ) 

- Q^{u x (u 3x + u x2y ) + u y {u 2xy + u 3y ))J +Vy^3(3u 2 u xy + 6^u xy (u 2x + u 2y )J = 0, (10) 
A^ 2 -fxM)+P x (t(|aM 3 -/3(^-^)+2/3M(M 2x + M 2y ))-|x(iaM 2 + /3M 2x ) + 2f^ 

+ Vy(- 2(3(tU X Uy+±XU X y)^) = 0. (11) 

Note that the fourth conservation law (iTTj) explicitly depends on t and x. 

Example 2. The well-known (2+l)-dimensional Kadomtsev-Petviashvili (KP) equation, 

(u t + auu x + u 3x ) x + a 2 u 2y = 0, (12) 

for u(x,y,t), describes shallow water waves with wavelengths much greater than their am- 
plitude moving in the x-direction a nd subject to weak variations in the y-direction 
( Kadomtsev and Petviashvilil . Il970l ). The parameter a occurs after a re-scaling of the phys- 
ical coefficients and a 2 = ±1. Obviously, the KP equation is not an evolution equation. 
However, it can be written as an evolution system in space variable y, 

u y = v, Vy = -cr 2 (u tx + au 2 x + auu 2x + u Ax ). (13) 

Note that -^ = a 2 , and thus cr 4 = 1. System (fl3"|) instead of (fl2"|) will be used in Section HI 
ConservationLawsMD.m has an algorithm that will identify an evolution variable and 
transform the given PDE into a system of evolution equations. 
Equation (I12p expresses conservation of momentum: 

V t {u x ) + V x (auu x + u 3x ) + Vy(a 2 Uy) = 0. (14) 

Other well- documented conservation laws ( WolJ . 120021 ) are 

Vtyfu) +V x (J(\au 2 + u 2x ) + {\a 2 f'y 2 - fx)(u t + auu x + u 3x )j 

+V y ((i/V - o 2 jx)u y - f'yu) = 0, (15) 

Vt\fyu\+V x (fy(\au 2 + u 2x ) + y{\a 2 fy 2 - fx)(u t + auu x + u 3x )j 

+V y (y(\fy 2 - a 2 fx)u y - (\fy 2 - a 2 fx)u^ = 0, (16) 

where / = f(t) is an arbitrary function. Thus, there is an infinite family of conservation 
laws, each of the form ffT5|) or ffT6|) . In Section H] we will show how fl8l)- ffTTl) are computed 
straightforwardly with our algorithm. We will also compute several conservation laws for the 
KP equation. Our current code does not (algorithmically) compute (IT51) and ffTBT) . Instead, 
conservation laws obtained with the code allow the user to conjecture and test the form of 
(fT5|) and ffTBT) . In Section [51 we give computational details and show how ffloT) and (fT5"j) can 
be verified. 
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3. Tools from the Calculus of Variations and Differential Geometry 

Three operators from the calculus of variations and differential geometry play a major 
role in the conservation law algorithm. Namely, the total derivative operator, and the Euler 
and homotopy operators. All three operators (which act on the jet space) can be defined 
algorithmically which allows for straightforward and efficient computations. 

The algorithm in Section H] requires that operations applied to differential functions 
/(x, u^ M ' ) (x)), take place in the jet space, where one component of x is a parameter. 

Although in later sections, one of the space variables will serve as the parameter, in this 
section we arbitrarily choose t as the parameter (matching (J5J)). Thus, in all definitions and 
theorems in this section, ID means that there is only one space variable, yet x = (x,t). 
Likewise, in 2 D and 3 D cases, x = (x, y, t) and x = (x, y, z, t), respectively. 

Using 05]), we assume that all partial derivatives of u with respect to t are eliminated 
from /. Thus, /(x,u^ M ^(x)) with 

u (x) = (u ,« I ,n !/ ,M 2 ,tt 2l ,?i 2!/ ,?i2 Z ,M ;E!/ ,...,« M jv lM w yM jv 2 ), (17) 

with M(, Ml-, Mg, and M as defined earlier. Each term in / must be a monomial in jet space 
variables, either multiplied with a constant or variable coefficient. 

Definition 1. The total derivative operator T> x in 2D is defined as 

df N Mi M * df 

j=lk 1 =0k 2 =0 ua k 1 xk 2 y 

where M( and M| are the orders of / for component m j with respect to x and y, respectively. 
T> y is defined analogously. Since t is parameter, T> t (in 2 D) is defined in a simpler manner, 

»*/=f+EEE^^>i (w) 

j=l k 1= k 2 =0 kixk 2 y 

If a total derivative operator were applied by hand to a differential function, /(x, u^ M ^(x)), 
one would use the product and chain rules to complete the computation. However, formulas 
like fll]), (ITB1) . and (IT^|) are more suitable for symbolic computation. 

The Euler operator (also kno wn as t he va riational derivative) plays a fundamental 



role in the calculus of variations (lOlverl . Il993l ). and serves as a key tool in our conser- 



vation laws algorithm. The Euler operator can be defined for any number of indepen- 
dent and dependent variables. For example in ID, the Euler operator is denoted by 

*-"u(a;) \'~'u 1 (x)j *-"u 2 (x)) • • • j *-"uJ(a:)j • • • i ^~'u N (x))- 



Definition 2. The ID Euler operator for dependent variable u J (x) is defined as 



fc=0 

df 



;_df_ 

du{ 



kx 



dip 



^ df - df „ df 



did 



du 3 



2x 



X du 3 3x 



■ ■ + {-V x ) M i 



df 



du 

M[x 



(20) 



j = 1 , . . . , N. The 2 D and 3 D Euler operators are defined analogously (jOlverl . 119931 ) . For 
example, the 2 D Euler operator is 



M{ M 3 2 

c uHx , y) f = ^2Y^(-v x ) k ^-v y ) k 

fc 1= 0fc 2 =0 



fc2 df 



du kl 



x k^y 



j = l,...,N. 



(21) 



The Euler operator allows one to test if differential functions are exact which is a key step 
in the computation of conservation laws. 

Definition 3. Let /be a differential function of order M. In 1 D, / is called exact if / is a 
total derivative, i.e., there exists a differential function -F(x, u^ M_1 ^(x)) such that / = T> X F. 
In 2D or 3D, / is exact if / is a total divergence, i.e., there exists a differential vector 
function F(x, u( M ~ 1 )(x)) such that / = DivF. 

Theorem 1. A differential function f is exact if and only if >C u (x)/ — 0. Here, is the 
vector (0, 0, ■ • • , 0) which has N components matching the number of components ofu. 



Proof. The proof for a general mult i- dimensional case is given in, e.g., iPoold (120091 ) 



Next, we turn to the homotopy operator (jAndersonl . l2004al ; lOlverl . 119931 ). which integrates 
exact 1 D differential functions, or inverts the total divergence of exact 2 D or 3 D differential 
functions. Integration routines in CAS have been unreliable when integrating exact differ- 
ential expressions involving unspecified functions. Often the built-in integration by parts 
routines fail when arbitrary functions appear in the integrand. The 1 D homotopy operator 
offers an attractive alternative since it circumvents integration by parts altogether. 



Definition 4. Let / be an exac t 1 D differential function. The homotopy operator in 1 D 
is defined (IHereman et all 120071 ) as 



K(x)/ = / I ^2l u i(x)f I [Au] 



d\ 



(22) 



where u = (u 1 , . . . ,vP, . . . ,u N ). The integrand, X u j^f, is defined as 



m{ 



'fc-i 



Zm(x)f = Y1 ^ u lx 



-V, 



vfe-(i+l) 



fc=l \j=0 






(23) 



where M\ is the order of / in dependent variable v? with respect to x. The notation /[Au] 
means that in / one replaces u by A u, u x by A u x , and so on for all derivatives of u. A is an 
auxiliary parameter that traces the homotopic path. 



Given an exact differential function, the 1 D homotopy operator f[2"2"j) replaces integration 
by parts (in x) with a sequence of differentiations followed by a standard integration with 
respect to A. Indeed, the following theorem states one purpose of the homotopy operator. 

Theorem 2. Let f be exact, i.e., T> X F = f for some differential function F(x, u^ M_1 ^(x)). 
Then,F = V- 1 f = n u{x) f. 



Proof. A proof for the 1 D cas e in the lang uage of standard calculus is given in 

Poole and Heremanl (2010J). See lOlverl f 19931 ) for a proof based on the variational complex. 



The homotopy operator (j22"l) has been a re liable tool f or int e grating exact polynomial dif- 
ferent ia l expressions. F o r applications, s ee Cheviakovi ( 20071 . |2010| ); iDeconinck and Nivala 
(120091 ); iHeremanl (120061 ); iHereman et all (120071 ) . However, the homot opy operator fails to 
i ntegr ate certain classes of exact rational expressi ons as discussed in iPoole and Hereman 
(2010). Although, the homotopy integrator code in IPoole and Heremanl ( 20091 ) covers large 
classes of exact rational functions, we will not consider rational expressions in this paper. 

CAS often cannot invert the divergences of exact 2 D and 3 D differential functions, 
although some capabilities exist in Maple. Again, the homotopy operator is a valuable tool 
to compute Div -1 , when it is impossible to do so by hand or by using the available software 
tools. 



Definition 5. The 2D homotopy operator is a "vector" operator with two components, 

where 
<»,/ = / (t 41,)/) M £ and H%J = jf (£, 41,,/) M y ■ ^ 



The x-integrand, 2^L \f, is given by 

M l M 2 /fci-l k 2 

'%J = E E £ E B< " <-» (-v,t-*-\-v y , 

fe 1= l fc 2 =0 \u=0 i 2 =0 
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k2-h 



Of 



du{ 



(26) 



k\xk,2y 



with combinatorial coefficient B^ x ' = B(ii,i 2 ,ki, k 2 ), where 

(h+i2\ (k 1 +k 2 -ii-i2~l\ 
\ ii ) \ k\—i\ — \ I 



dcf 

Similarly, the y- integrand, 1^, -,/, is defined as 



b^ 2 mmv= — '\ kl xr • (27) 



: (s/) 

W (a 
M l M 2 / fci fc 2 -l 



<!,)/ = E E E E ^ «Uy (-V x ) k ^(-V y f^ -?—, (28) 

fc 1= 0fc 2 =l \n=0i 2 =0 / OU k 1 xk 2 y 

where B^> = B(i 2 , ii, k 2 , k\). 



Definition 6. The homotopy operator in 3 D is a three-component vector operator, 

\ K ™u(x,y,z)J>™u(x,y,z)J>™-u(x,y,z)-J J ' (^9) 

where the x-component is given by 



^Cow)/ - / I z2 X ^(x, y ,z)f J t Au l 



dX 

T' 



(30) 



The y- and ^-components are defined analogously. The x-integrand is given by 

M{ M| M{ kl -l k 2 k 3 

ui(x,y,z)J = 2-^ Z_^ 2-^ 2-^ 2-^2-^t V U hxi 2 yi 3 z 

fc 1 = l k 2 =0 k 3 =0 h=0 i 2 =0 i 3 =0 

{ _ Vx) k.-n-l { - Vy f^ {-V z f^) — ?l , (31) 

k\xk 2 y k 3 z 

with combinatorial coefficient B^ x > = B(ii,i 2 ,i3,ki,k 2 , k 3 ) where 

li\+i 2 +i 3 \ li 2 +i 3 \ (ki+k 2 +k 3 -i 1 -i 2 -i 3 -l\ lk 2 +k 3 -i 2 -i 3 \ 

^^1^2,^3,^1,^2,%; /k 1+ k 2 +k 3 \ (k 2 +k 3 \ • \ 6Z > 



The integrands I 2 , f and I /, , f are defined analogously. Based on cyclic per- 

' u3(x,y,z) J u3 (x,y,z) J ■■ ■ ' 



r <i J and M 

u3(x,y,z) J ui(i^ 

mutations, they have combinatorial coefficients B^ y > = B(i 2 , 13, i\, k 2 , k^, k\) and B^' 
B(i 3 , ii, i 2 , k 3 , fci, k 2 ), respectively. 

Using homotopy operators, Div -1 can be computed based on the following theorem. 

Theorem 3. Let f be exact, i.e., f = DivF for some F(x, ir M_1 )(x)). Then, in the 2D 
case, F = Div^ 1 / = \U^ y) f,U^ y) fj ■ Analogously, in 3D one has 

f = Div- 1 / = (Vf j,n {y ) j,n (z ) j). 

■1 \ u(x,y,z) J ' u(x,y,z) J ' u(x,y,z) J I 



Proof. A proof for the 2D case is given in iPoold (120091 ). The 3D case could be proven 
with similar arguments. 

Unfortunately, the outcome of the homotopy operator is not unique. The homotopy in- 
tegral in the ID case has a harmless arbitrary constant. However, in the 2D and 3D 
cases there are infinitely many non-trivial choices for F. From vector calculus we know 
that DivCurlK = 0. Thus, the addition of CurlK to F would not alter DivF. More pre- 
cisely, for K = (V y 9, -vj) in 2 D, or for K = (V y n - V z £, V z 9 - V x rj, V x £ - V y 9) in 3 D, 
Div G = Div (F + K ) = DivF, where 9, r j, and £ are arbitrary functions. To obtain a con- 
cise result for Div -1 , IPoole and Hereman ( 20101) developed an algo rithm that removes curl 



terms. Furthermore, when / is rational ( IPoole and Heremanl . 120101 ). the homotopy operator 



may fail at the singularities of /; but rational functions are not considered in this paper. 

4. An Algorithm for Computing a Conservation Law 

To compute a conservation law, the PDE is assumed to be in the form given in (JSJ) for 
a suitable evolution variable. Adhering to ([3]), if the evolution variable is t, we construct a 
candidate density. However, if the evolution variable is x, y, or z, we construct a candidate 
component of the flux corresponding to the evolution variable. For argument's sake let us 
assume that the evolution variable is time. 

The candidate density is constructed by taking a linear combination (with undetermined 
coefficients) of terms that are invariant under the scaling symmetry of the PDE. The total 
time derivative of the candidate is computed and evaluated on fl5]), thus removing all time 
derivatives from the problem. The resulting expression must be exact, so we use the Euler 
operator and Theorem[T]to derive the linear system that yields the undetermined coefficients. 
Substituting these coefficients into the candidate leads to a valid density. 

Once the density is known the homotopy operator and Theorems [2] or [3] are used to 
compute the associated flux, J, taking advantage of (J3J). 

In contrast to other algorithms which attempt to compute the components of P in (J2]) all 
at once, our algorithm computes the density first, followed by the flux. Although restricted 
to polynomial conservation laws, our constructive method leads to short densities (which 
are free of divergences and divergence-equivalent terms) and fluxes in which all curl terms 
are automatically removed. 

Definition 7. A term or expression / is a divergence if there exists a vector F such that 
/ = DivF. In the ID case, / is a total derivative if there exists a function F such that 
/ = T> X F. Note that T> x f is essentially a one-dimensional divergence. So, from here onwards, 
the term "divergence" will also cover the "total derivative" case. Two or more terms are 
divergence-equivalent when a linear combination of the terms is a divergence. 

To illustrate the subtleties of the algorithm we intersperse the steps of the algorithm with 
two examples, viz., the ZK and KP equations. 
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4-1. Computing the Scaling Symmetry 

A PDE has a unique set of Lie-point symmetries whi ch may include trans lations, rota- 
tions, dilations, Galilean boosts, and other symmetries ( iBfuman et all l2010l ). The appli- 
cation of such symmetries allows one to generate new solutions from known solutions. We 
will use only one type of Lie-point symmetry, namely, the scaling or dilation symmetry, to 
formulate a "candidate density." 

Let us assume that a PDE has a scaling symmetry. For example, the ZK equation (j7|) is 
invariant under the scaling symmetry 

(x, y, t, u) ->■ (X~ x x, X^y, \~H, X 2 u), (33) 

where A is an arbitrary scaling parameter, not to be confused with A in Definitions 0] 
through [BJ 

Step 1-ZK (Computing the scaling symmetry). To compute (133]) with linear algebra, 
assume that (J7|) for u(x, y, t) scales uniformly under 

(x, y, t, u) -+ (X, Y, T, U) = (X a x, X b y, A% X d u), (34) 

where U(X, Y, T) and a, b, c, and d are undetermined (rational) exponents. We assume that 
the parameters a and j3 do not scale. By the chain rule, ([7]) transforms into 

u t + auu x + f3u 3x + (3u x2y 

= X c - d (U T + aX a - c - d UU x + fiX 3a - c U 3X + PX a+2b - c U X2Y ) = 0. (35) 

If a — c — d = 3a — c = a + 2b — c = 0, we have (J7J) for U(X, Y, T) up to the scaling factor 
X c ~ d . Setting a = —1, we find b = —1, c = —3, and d = 2, corresponding to (1331 . 

Step 1-KP (Computing the scaling symmetry). The scaling symmetry for the KP equation 
will be computed similarly. Assume that ( fl~3]) scales uniformly under 

(x, y, t, u, v) -> (X, Y, T, U, V) = (X a x, X b y, A c t, X d u, X e v) , (36) 



with unknown rational exponents a through e. Applying the chain rule to get (113]) expressed 
in the variables (X, Y, T, U, V) yields 

Uy - v = X h - d {U Y - X d ~ b - e V) = 0, 
v y + a 2 {u tx + ul + uu 2x + u ix ) 

= X b ~ e (V Y + a 2 (X a - b+c - d+e U TX + aX 2a ~ b ~ 2d+e (U x + UU 2X ) + X ia - b - d+e U iX )) = 0.(37) 



By setting d — b — e = a — b + c — d + e = 2a — b — 2d + e = Aa — b — d + e = 0, (1371) becomes a 
scaled version of (TT3]) in the new variables U(X, Y, T) and V(X, Y, T). Setting a = — 1 yields 
b = —2, c = —3, d = 2, and e = 4. Hence, 

(x, y, t, u, v) ->■ (A _1 x, A~ 2 y, A~ 3 t, X 2 u, AS) (38) 

is a scaling symmetry of (TT3]) . 
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4-2. Constructing a Candidate Component 

Conservation law (J2J) must hold on solutions of the PDE. Therefore, we search for poly- 
nomial conservation laws that obey the scaling symmetry of the PDE. Indeed, we have yet 
to find a polynomial conservation law that does not adhere to the scaling symmetry. 

Based on the scaling symmetry of the PDE, we choose a scaling factor for one of the 
components of P in 02]). The selected scaling factor will be called the rank (R) of that 
component. Then, we construct a candidate for that component as a linear combination of 
monomial terms (all of rank R) with undetermined coefficients. By dynamically removing 
divergence terms and divergence-equivalent terms that candidate is short and of low order. 

Step 2-ZK (Building the candidate component). Since the ZK equation ([7|) has t as evolu- 
tion variable, we will compute the density p of ([3]) of a fixed rank, for example, R = 6. 

(a) Construct a list, V, of differential terms containing all powers of dependent variables 
and products of dependent variables that have rank 6 or less. By (13"3"]) . u has a scaling factor 
of 2, so u 3 scales to rank 6 and u 2 has rank 4. This leads to V = {u 3 , u 2 , u}. 

(b) Bring all of the terms in V up to rank 6 and put them into a new list, Q. This is done 
by applying the total derivative operators with respect to the space variables. Taking the 
terms in V, u 3 has rank 6 and is placed directly into Q. The term u 2 has rank 4 and can be 
brought up to rank six in three ways: either by applying D x twice, by applying T> y twice, 
or by applying each of T> x and T> y once, since both T> x andP y have scaling factors of 1. All 
three possibilities are considered and the resulting terms are put into Q. Similarly, the term 
u can be brought up to rank 6 in five ways, and all results are placed into Q. Doing so, 

Q = {u 3 , u 2 x , uu 2x , u 2 y , uu 2y , u x u y , uu xy , u Ax , u 3xy , u 2x2y , u x3y , u Ay ], (39) 

in which all monomials are now of rank 6. 

(c) With the goal of constructing a nontrivial density with the least number of terms, remove 
all terms that are divergences or are divergence- equivalent to other terms in Q. This can be 
done algorithmically by applying the Euler operator ( 12 ip to each term in ( I39p . yielding 

Cu(x, y )Q = {3u 2 ,-2u 2x ,2u 2x , -2u 2 y,2u 2y ,-2u X y,2u xy , 0,0, 0,0,0}. (40) 

By Theorem [U divergences are terms corresponding to in ( 14"U|) . Hence, u Ax , u 3xy , u 2x2y , 
u x3y , and u Ay are divergences and can be removed from Q. Next, all divergence-equivalent 



terms will be removed. Following iHereman et al\ ( 120051 ). form a linear combination of the 



terms that remained in (140]) with undetermined coefficients pi, gather like terms, and set it 
identically equal to zero, 

3piu 2 + 2(p 3 - p 2 )u 2x + 2(p 5 - p 4 )u 2 y + 2(p 7 - p 6 )u xy = 0. (41) 

Hence, p\ = 0, p 2 = p 3 , p A = p 5 , and p 6 = p 7 . Terms with coefficients p 3 , p 5 , and p 7 are 
divergence-equivalent to the terms with coefficients p 2 . p A , and p 6 , respectively. For each 
divergence-equivalent pair, the terms of highest order are removed from Q in (I3"9"]) . After all 
divergences and divergence-equivalent terms are removed, Q = {u 3 , u^u 2 , u x u y }. 
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(d) A candidate density is obtained by forming a linear combination of the remaining terms 
in Q using undetermined coefficients q. Thus, the candidate density of rank 6 for ([7]) is 



P 



Ci"U 3 + C 2 u\ + CsUy + c^u x u y . 



(42) 



Now, we turn to the KP equation (fT2"j) . The conservation laws for the KP equation, 
(|15j) and (fTol) . involve an arbitrary functional coefficient /(£). The scaling factor for f(t) 
depends on the degree if f(t) is polynomial; whereas there is no scaling factor if f(t) is 
non-polynomial. In general, working with undetermined functional (instea d of c o nstan t) 
coefficients f(x, y, z, t) would require a sophisticated solver for PDEs for / (see lWolJ (120021 ) ). 
Therefore, we can not automatically compute (IT5|) and (II 6p with our method. However, our 
algorithm can find conservation laws with explicit variable coefficients, e.g., tx 2 ,txy, etc., as 
long as the degree is specified. Allowing such coefficients causes the candidate component 
to have a negative rank. By computing several conservation laws with explicit variable 
coefficients it is possible (by pattern matching) to guess and subsequently test the form of 
a conservation law with arbitrary functional coefficients. 

Step 2-KP (Building a candidate y-component). When the KP equation is replaced by 
(113)1 . the evolution variable is y. Thus, we will compute a candidate for the y-component of 
the flux, J y , in ([3]). The y-component will have rank equal to —3. The negative rank occurs 
since differential terms for the component are multiplied by Cit m x n y p , which, by (138]) . scales 
with A -3m A~ n A~ 2p = \-( 3m + n + 2 p) 5 where m, n, and p are positive integers. The total degree 
of the variable coefficient t m x n y p , is restricted to0<m + n + p<3. 

(a) As shown in Table [TJ construct two lists, one with all possible coefficients t m x n y p up to 
degree 3 and the other with differential terms, organized so that the combined rank equals 
—3. The rank of each term is computed using the scaling factors from ( 138]) . For example, t 
and x have scaling factors of —3 and — 1, respectively, so tx 2 has rank —5. Variable u has 
scaling factor 2, so t 2 xu has rank —3. Since we are computing the y-component of J, the 
differential terms contain only derivatives with respect to x and t. 



Factors of Type t m x n y p 


Differential Terms 


Product 


Rank 


Coefficient 


Rank 


Term 


Rank 


-5 


tx 2 , xy 2 , ty 


2 


u 


-3 


-6 


y 3 , txy, t 2 


3 


Ux 


-3 


-7 


t 2 x, ty 2 


4 


u 2 , U 2x , V 


-3 


-8 


t 2 y 


5 


uu x , u t , u 3x , v x 


-3 


-9 


t 3 


6 


M 3 , UV, U 2 X , Utx, UU 2x , U 4x , V 2x 


-3 



Table 1: 
-3. 



Factors t m x n y p of degree 3 are paired with differential terms so that their products have ranks 
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(b) Combine the terms in Tabled] to create a list of all possible terms with rank —3, 

Q = {tx 2 u, xy 2 u, tyu, y 3 u Xl txyu x , t 2 u x , t 2 xu 2 ,ty 2 u 2 , t 2 xu 2x , ty 2 u 2x , t 2 xv, ty 2 v 

t 2 yuu x , t 2 yu t , t 2 yu 3x , t 2 yv x , t 3 u 3 , t 3 uv, t 3 u 2 x , t 3 u tx , t 3 uu 2x , t 3 u ix , t 3 v 2x }. (43) 

(c) Remove all divergences and divergence-equivalent terms. Apply the Euler operator to 
each term in (j43p . Next, linearly combine the resulting terms to get 

tx 2 \ fxy 2 \ fty\ fty\ (2t 2 xu\ (2ty 2 u\ ( 



, /0\ f2ty\, /3t 3 u 2 \ ft 3 v\ (2t 3 u 2x \^ (2t 3 u 2x \ 

+Pl2 [ty 2 )- pu [o ) +Pl7 { ) +P ^[t 3 u)' Pw { ) +p *{ J=M 44 ) 

where the subscript of the undetermined coefficient, Pi, corresponds to the ith term in Q. 
Missing Pi correspond to terms that are divergences. Gather like terms, set their coefficients 
equal to zero, and solve the resulting linear system for the p^, to get p\ = p 2 = p 7 = p 8 = 
Pn = P12 — Pn — Pis — 0, P3 = Ps + 2pi4, and pi$ = p 2 \. Thus, both terms with coefficients 
P5 and P14 are divergence-equivalent to the term with coefficient p 3 . Likewise, the term with 
coefficient p 2 \ is divergence-equivalent to the term with coefficient pig. For each divergence- 
equivalent pair, the terms with the highest order are removed from (|4~3]) . After removal of 
divergences and divergence-equivalent terms 

Q = {tx 2 u, xy 2 u, tyu, t 2 xu 2 , ty 2 u 2 , t 2 xv, ty 2 v, t 3 u 3 , t 3 uv, t 3 u 2 x }. (45) 

(d) A linear combination of the terms in (145 j) with undetermined coefficients q yields the 
candidate (of rank —3) for the ^-component of the flux, i.e., 

J y = C\tx 2 u + c 2 xy 2 u + c 3 tyu + c^xu 2 + c 5 ty 2 u 2 + c 6 t 2 xv 

+ c 7 ty 2 v + c 8 t 3 u 3 + c Q t 3 uv + c w t 3 u 2 x . (46) 

4-3. Evaluating the Undetermined Coefficients 

All, part, or none of the candidate density fj42|) may be an actual density for the ZK 
equation. It is also possible that the candidate is a linear combination of two or more inde- 
pendent densities, yielding independent conservation laws. The true nature of the density 
will be revealed by computing the undetermined coefficients. By 03]), T> t p = — Div (J x , J y ), 
so T> t p must be a divergence with respect to the space variables x and y. Using Theorem [1] 
an algorithm for computing the undetermined coefficients readily follows. 

Step 3-ZK (Computing the undetermined coefficients). To compute the undetermined 
coefficients, we form a system of linear equations for these coefficients. As part of the 
solution process, we also generate compatibility conditions for the constant parameters in 
the PDE, if present. 

(a) Compute the total derivative with respect to t of (j4"2"|) . 

V t p = 3ciu 2 u t + 2c 2 u x u tx + 2c 3 u y uty + c^u tx u y + u x u ty ). (47) 
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Let E = —V t p after u t and u tx have been replaced using (J7|). This yields 

E = 3ciu 2 (auu x + /3(u 3x + u x2y )) + 2c 2 u x (auu x + f3(u 3x + u x2y )) x 

+ 2c 3 u y (auu x + f3(u 3x + u x2y )) y + c 4 (u y (auu x + (3(u 3x + u x2y )) x 

+ u x (auu x + (3(u 3x + u x2y ))y). (48) 

(b) By 0, E = Div (J x ,jy). Therefore, by Theorem QJ C u {x, y ) E = 0. Apply the Euler 
operator to (148]) . gather like terms, and set the result identically equal to zero: 



= Cu(x, v ) E = -2((3ci/3 + c 3 a)u x u 2y + 2(3c 1 (3 + c 3 a)u y 



u 



■>■(/ 



+ 2c 4 au x u xy + c 4 au y u 2x + 3(3ci/3 + c 2 a)u x u 2x ) . (49) 

(c) Form a linear system of equations for the undetermined coefficients q by setting each 
coefficient equal to zero, thus satisfying (f49|) . After eliminating duplicate equations, the 
system is 

3cx/? + c 3 a = 0, c 4 a = 0, 3ci/3 + c 2 a = 0. (50) 

(d) Check for possible compatibility conditions on the parameters a and j3 in floUl) . This 
is done by setting each q = 1, one at a time, and alg e braic ally eliminating the other unde- 



termined coefficients. Consult iGoktas and Heremanl (119971 ) for details about searching for 



compatibility conditions. System (150]) is compatible for all nonzero a and j3. 

(e) Solve (150]) . taking into account the compatibility conditions (if applicable). Here, 

c 2 = c 3 = -3|ci, c 4 = 0, (51) 

where c\ is arbitrary. We set C\ — 1 so that the density is normalized on the highest degree 
term, yielding 

P = u'-^M + ul). (52) 

Step 3-KP (Computing the undetermined coefficients). The procedure to find the unde- 
termined coefficients in the KP case is similar to that of the ZK case. 

(a) Starting from (T46]) . compute 

T> y J y = (c\x + 2c4tu)txu y + c 2 xy(2u + yu y ) + (c 3 t + 2c^,tyu){u + yu y ) 

+ c e t 2 xv y + c 7 ty{2v + yVy) + 3c 8 t 3 u 2 u y + c 9 t 3 (u y v + uv y ) + 2c w t 3 u x u xy , (53) 

and replace u y and v y and their differential consequences using (ITBl . Thus, 

E = —V y J v = —{c\X + 2citu)txv — c 2 xy(2u + yv) — (c 3 t + 2c 5 tyu)(u + yv) 

+ a (cet x + city + cgt u)(u tx + au x + ctuu 2x + u ix ) — 2cjtyv 

- 3c 8 t 3 u 2 v - c 9 t 3 v 2 - 2ci t 3 u x v x . (54) 
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(b) Apply the Euler operator to (1511) and set the result identically equal to zero. This yields 

(0, 0) = 0=£ u ( t)a ;) E— i£ u (t jX ) E, C v (t, x ) EJ 

= — ( 2c 2 xy + (c 3 — 2cr 2 c 6 )t + 2c 4 t 2 xv + 2c 5 ty(2u + yv) + Qc s t 2 uv 

— 2cr 2 c 9 £ 2 (|a E + tu tx + atu 2 x + auu 2x + tu± x ) — 2c w t 3 v 2x , C\tx 2 + c 2 xy 2 

+ (c 3 + 2c 7 )ty + 2c 4 t 2 xu + 2c 5 ty 2 u + 3c 8 t 3 u 2 + 2c 9 t 3 v - 2ci t 3 u 2 < x \ ■ (55) 

(c) Form a linear system for the undetermined coefficients q. After duplicate equations and 
common factors have been removed, one gets 

Ci = 0, c 2 = 0, c 3 - 2a 2 c 6 = 0, c 3 + 2c 7 = 0, c 4 = 0, c 5 = 0, c 8 = 0, c 9 = 0, c 10 = 0. (56) 

(d) Compute potential compatibility conditions on the parameters a and a. Again, the 
system is compatible for all nonzero values of a and a. 

(e) Use a 2 = ±1 and solve the linear system, yielding 

ci = c 2 = c 4 = c 5 = c 8 = c 9 = cio = 0, c 6 = 7jO- 2 c 3 , c 7 = -|c 3 . (57) 

Set c 3 = —2 (to normalize the density) and substitute the result into (146 p . to obtain 

jy = -t(2yu + (crHx - y 2 )v), (58) 

which matches J y in ffi~5]) if f(t) = t 2 and v = u y . 

4-4- Completing the Conservation Law 

With the density (or a component of the flux at hand), the remaining components of the 
conservation law can be computed with the homotopy operator using Theorem [5] or [3J 

Step 4-ZK (Computing the flux, J). Again, by the continuity equation (j^J), DivJ = 
Div ( J x , J v ) = —T> t p = E. Therefore, compute Div -1 E, where the divergence is with respect 
to x and y. After substitution of ( 15 ip with c\ — \ into 



E = 3u 2 (auu x + j3u Zx + /3u x2y ) - 6^u x (auu x + (3u 3x + f3u x2y ) x 

- 6^u y (auu x + (3u 3x + /3u x2y ) y . (59) 

Apply the 2 D homotopy operator from Theorem [3j Compute the integrands (126]) and (1281) : 



/g iy) E = 3au 4 + $ {9u 2 (u 2x + \u 2v ) - 6u(3« 2 + u 2 y )) + f {6u 2 2x + 5u 2 xy + \u 2 2y 

+ \u(U 2x2 y + U A y) - U X (12U 3X + 7^2^) - Uy(3U 3 y + 8U 2x y) + ^U 2x U 2 y) , (60) 

Mx.j/)^ = 3(3u(uu xy - Au x u y ) - |^- (3m(m 3:C j / + w^) + M a .(13M 2a: y + 3w % ) 

+ 5m^(m 3x . + 3^22/) - 9«zj,(w2s + W 2 j/)) , (61) 

16 



respectively. Use fl25|) . to compute J = \H*L V )E,'H^L V \E\ where 

J 



dX 

T 
f mz 4 + {3u 2 (u 2x + |u 22/ ) - 2u(3u 2 x + u 2 y )) + £ (3< + f<, + |^ 
+ \u(u 2x2y + u 4y ) - u x (6u 3x + \u x2y ) - u y (\u 3y + Au 2xy ) + \u 2x u 2y ) , (62) 



n {y ] ,E = /Yx? ^ [Au] ^ 

u(a;,2/) / ^ u(a;,y) y L J \ 



= (3u(uu xy - Au x u y ) - \^ {3u(u 3xy + u x3y ) + u x (13u 2xy + 3u 3y ) 

+ 5u y (u 3x + 3u x2y ) - 9u xy (u 2x + u 2y )) . (63) 

Notice that J has a curl term, K = (T> y 8, —D x 9), with 

9 = 2/3u 2 u y + \^ [3u(u 2xy + u 3y ) + 5(2u x u xy + 3u y u 2y + u 2x u y ) ) . (64) 

Therefore, compute J — K to obtain 

,F = 3(u 2 {\au 2 + (3u 2x ) - 2(3u(u 2 x + u 2 y ) + ^[u\ x - u\ y ) 

- 2^(u x (u 3x + u x2y ) + u y (u 2xy + u 3y ))\ , (65) 

jy = 3/3 (u 2 u xy + 2^u xy (u 2x + u 2y )\ , (66) 

which match the components in (1101) . 

Step 4-KP (Computing the density and the x-component of the flux). For the KP example, 
(p, J x ) remains to be computed. Using the continuity equation ([3]), V t p + V x J x = —V y J y = 
E. Thus, to find (p, J x ), compute Div~ 1 E, where this time the divergence is with respect to 
t and x. Proceed as in the previous example. First, substitute (15TJ) and c 3 = —2 into ( 15"4"|) . 

E = t(2u+ (a 2 y 2 - tx) (u tx + au 2 x + auu 2x + u ix ) J . (67) 

Second, compute the integrands for the homotopy operator, 

l{ X x) E = -\(W X - u x Z) — = \t{tu + [a 2 y 2 - tx)u x ), (68) 

I%, X) E = 0, (69) 

t(x) J-, 9E , , . oE . oE . n o ^r-, t-\ uE 
tutt x) E = M 7i 2 ( uT> t - u t X ) 7^ \ uV x - U * X ) 7^ \ uV x - u xV x + u 2x V x - u 3x X} 



lift <r1 <~1 o\ u,J -^t U 't- Lj I r\ V S '-"X-^ I r\ \ X "I I TU, »I I Uj 3 X-L-J o 

"^ du x 2 du tx du 2x du± x 

= t 2 (au 2 + u 2x ) + t(a 2 y 2 - tx)(\u t + 2aMM :I . + -u 3x ) - (|cr 2 y 2 - tx)u, (70) 

4t)^ = °- ( 71 ) 
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Next, compute 

P = K% >X) E = f (lg x) E + lil x) E) [Xu] 

(|t(tu + [p 2 y 2 — tx)u x )^j dX = ^t(tu + (a 2 y 2 — tx)u x ), 

/ (t 2 (aXu 2 + u 2x ) + t(a 2 y 2 — tx){\u t + 2aXuu x + u 3x ) — 
Jo 



j x = -h { :L) e 



,dX 

T 



(72) 



V- 2 



y 2 — tx)u) dX 



= t 2 (\au 2 + u 2x ) + t(a 2 y 2 - tx)(\u t + auu x + u 3x ) - i\o 2 y 2 - tx)u, 
and remove the curl term K = (V x 9, —V t 6) with 9 = \t{o 2 y 2 — tx)u, to obtain 



(73) 



(74) 



p = t u, J x = t (^au + u 2x ) + t(o y — tx)(u t + auu x + u 3x ). 
The computed conservation law is the same as (JT5J) where f(t) = t 2 and v = u y . 

5. A Generalized Conservation Law for the KP Equation 

Due to the presence of an arbitrary function f(t), it is impossible to algorithmically 
compute f TTo^) with our code. The generalization of (|7I|) to ([To'j) is based on inspection of 
the conservation laws in Table |2] as computed by our program ConservationLawsMD.m. 
Indeed, pattern matching with the results in Tableland some interactive work lead to (TT5l) . 



Rank 


Conservation Law 


5 


vAuj + V X ( \au 2 + u 2x - x(u t + auu x + u 3x ) j + V y ( - a 2 xv) = 


2 


V t \tu) + V x U(\au 2 + u 2x ) + {\o 2 y 2 - tx){u t + auu x + u 3x )\ 
+ V y ({\y 2 - a 2 tx)v - yu) = 0. 


-1 


V t \t 2 u\ + V x (t 2 (\au 2 + u 2x ) + t(a 2 y 2 - tx)(u t + auu x + u 3x )\ 
+ V y (t(y 2 - a 2 tx)v - 2tyu) = 


-4 


V t (t 3 u) + £> x .(V(|cm 2 + u 2x ) + t 2 (\a 2 y 2 - tx)(u t + auu x + u 3x )) 
+ Vy(t 2 ((^y 2 -a 2 tx)v-3yuj) =0 



Table 2: Additional conservation laws for the KP equation ([HI)) . 

which can be then be verified with ConservationLawsMD.m as follows. 

The conservation laws in Table [2] suggest that a density has the form t n u, or more general, 
f{t)u, where f(t) is an arbitrary function. The corresponding flux would be harder to guess. 
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However, it can be computed as follows. Since the KP equation ( TT3|) is an evolution equation 
in y, we construct a suitable candidate for J y . Guided by the results in Table [2j we take 

jy = c x f{t)yu + c 2 f(t)y 2 v + c 3 f(t)xv, (75) 

where Ci, C2, and c 3 are undetermined coefficients, and u y is replaced by v in agreement with 
(1T3|) . As before, we compute T> y J y and replace u y and v y using ( 113]) . Doing so, 

i? = Py-P 27 = ci/'w + (ci + 2c 2 )f'yv - {a 2 c 2 f'y 2 + a 2 c 3 fx)(u tx + au 2 x + auu 2x + u^)- (76) 

By ©, Vyjy = -Div(p, J x ). By Theorem [TJ 

(0, 0) = = £ u(t)Cc) £ = (( Cl - a 2 c 3 )/', ( Cl + 2c 2 )/V) • (77) 

Clearly, c 2 = —\c\ and c 3 = o 2 C\. If we set c\ = —1 and v = u y we obtain J y in ([TBI) . 
Application of the homotopy operator (in this case to an expression with arbitrary functional 
coefficients) yields (p,J x ). This is how conservation law (TT5l) was computed. Conservation 
law (|T6|) was obtained in a similar way. Both conservation laws were then verified using the 
ConservationLawsMD.m code. 

6. Applications 

In this section we state results obtained by using our algorithm on a variety of (2+1)- and 
(3+l)-dimensional nonlinear PDEs. The selected PDEs highlight several of the issues that 
arise when using our algorithm and software package ConservationLawsMD.m. 

6.1. The Sawada-Kotera Equation in 2D 



The (2+l)-dimensional SK equation (IKonopelchenko and Dubrovskyl . Il984l ). 

u t = 5u 2 u x + 5uu 3x + 5uu y + 5u x u 2x + 5u 2xy + u 5x - bd~ l u 2y + hUxd^Uy^ (7£ 
with u (x) = u(x, y, t) is a completely integrable 2 D generalization of t he standard SK equa- 



tion. The latter has infinitely many conservation laws (see, e.g., iGoktas and Hereman 



(Il997l )). Our algorithm can not handle the integral terms in (1751) . so we set v = d x 1 u 



x "V 



Doing so, (1751) becomes a system of evolution equations in y: 

v y = -~u t + u 2 u x + uu 3x + uv x + u x u 2x + v 3x + \u 5x + u x v, u y = v x . (79) 

Application of our algorithm to (1791) yields several conservation laws, all of which have 
densities 14, tu, t 2 u, etc., and yu, tyu, t 2 yu, etc. Like with the KP equation, this suggests 
that there are conservation laws with an arbitrary functional coefficient /(£). Proceeding as 
in Section [5] and using ConservationLawsMD.m, we obtained 

V t [fuj +V X (j'yv- 5/(|iz 3 + uv + uu 2x + u xy + \u Ax )) +V y (hfv - /V) = 0, (80) 

A (fyu) +V X \A\f'y 2 - 5fx)v- 5/y(|u 3 + uv + uu 2x + u xy + \u Ax )\ 

+V y (5fyv - {\fy 2 - 5fx)u) = 0. (81) 
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Note that the densities in ( 1HU1) and (IHTi) are identical to those in (|T5|) and (|T6|) for the 
KP equation. These two densities occur often in (2+l)-dimensional PDEs that have a Ut x 
instead of a Ut term, as shown in the next example. 

6.2. The Khokhlov-Zabolotskaya Equation in 2D and 3D 

The Khokhlov-Zabolotskaya (KZ) equation or dispersionless KP equ ation describes the 



propag ation of sound in non-linear media in two or three space dimensions (ISanders and Wang 



1997aj). The (2+l)-dimensional KZ equation, 

(u t - uu x ) x - u 2y = 0, (82) 

with m(x) = u(x,y,t) can be written as a system of evolution equations in y, 

u y = v , Vy = u tx -u 2 x - uu 2x , (83) 

by setting v = u y . Again, two familiar densities appear in the following conservation laws, 
computed indirectly as we showed for the KP and SK equations, 

V t (u x )+V x (-uu x )+V y (-u y ) = 0, (84) 

V t (fu)+V x (-\fu 2 -{\fy 2 + fx)(ut - uu x ))+V y ((\fy 2 + fx)u y -fyu) = 0, (85) 

V t (fyu)+V x (- \fyu 2 -y{\fy 2 + fx)(u t - uu x )) 

+V y (y(lfy 2 + f x ) Uy -(\fy 2 + fx)u) = 0, (86) 



where f(t) is an arbitrary function. Actually, (l8"o"|) and (I8"B"|) are nonlocal because, from 

u t — uu x = J u 2y dx. By swapping terms in the density and the x-component of the flux, 

( 185]) with f(t) = 1, can be rewritten as 



V t (xu x ) + V x Um 2 - xuu x ) +T> y (- xu y ) = 0, (87) 

which is local. The computation of conservation laws for the (3+l)-dimensional KZ equation, 

(u t - uu x ) x - u 2y - u 2z = 0, (88) 

where u(x) = u(x,y,z,t), is more difficult. This equation can be written as a system of 
evolution equations in either y or z. Although the intermediate results differ, either choice 
leads to equivalent conservation laws. Writing ( 18"B"|) as an evolution system in z, 

u z = v, v z = u tx -u 2 x - uu 2x - u 2y , (89) 

ConservationLawsMD.m is able to compute a variety of conservation laws whose densi- 
ties are shown in Table |3j 
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Rank 


Densities Explicitly Dependent on x, y, z 


2 


Pi = xu x 





92 = xyu x , p 3 = xzu x 


-1 


P4 = tU 


-2 


p 5 = xyzu x , p e = x{y 2 - z 2 )u x 


-3 


p 7 = tyu, p 8 = tzu 


-4 


p 9 = t 2 u, p w = xy(y 2 - 3z 2 )u x , p u = xz 2 (3y - z)u x 


-5 


P12 = tyzu, pi 3 = t(y 2 - z 2 )u 


-6 


Pu = t 2 yu, P15 = t 2 xzu, pi 6 = xyz(y 2 - z 2 )u x , p 17 = x(y i - 6y 2 z 2 


+ z 4 )u x 


-7 


pis = t 3 u, pig = ty(y 2 - 3z 2 )u, p 20 = tz(3y 2 - z 2 )u 


-8 


P2i = t 2 yzu, p 2 2 = t 2 (y 2 - z 2 )u, p 23 = xy(y 4 -10y 2 z 2 + 5z 4 )u x , 
p 2 4 = xz(5y 4 - 10y 2 z 2 + z 4 )u x 



Table 3: Densities for the (3+l)-dimensional KZ equation 



Density pi = xu x in Table [3] is part of local conservation law 



Vt xu x 



T> x [ \v 2 — xuu s 



X\ 






+ V 2 



xu. 



(90) 



t \ ^ ^x I i -^ x \ 9 w • XjlA ' U'x I i i ^ y \ ^ ^y 

which can be rewritten as a nonlocal conservation law 

V t (u) +Vj-\u 2 - x(u t - uu x )\ +Vy(xuyj +V z [xu z ) = 0. (91) 

In general, if a factor xu x appears in a density then that factor can be replaced by u. Doing 
so, all densities in Table[3]that can be expressed as p = g(y, z, t)u, where g(y, z, t) is arbitrary. 
Introducing an arbitrary function h = h(y, z, t), the conservation laws corresponding to the 
densities in Table [3] can be summarized as 

vAgvA+vJ- \gu 2 - (xg + h){u t - uu x )\+V y \Axg + h)u y - (xg y + h y )u 

+V Z ((xg + h)u z - (xg z + h z )uj = - [h 2y + h 2z ~g t + x{g 2y + g2z))u. (92) 

Equation (192]) is only a conservation law when the constraints Ag = and Ah = gt are 
satisfied, where A = -^ + ^. Thus, g must be a harmonic function and h must satisfy the 
Poisson equation with g t on the right hand side. Combining both equation s produces the 



eqi 
biharmonic equation A 2 h = 0. As shown by iTikhonov and Samarskiil (J1963I ). A 2 h = has 
general solutions of the form 

h = yh 1 (y,z) + h 2 (y,z) and h = zhi(y, z) + h 2 (y, z), (93) 

where Ahi = and Ah 2 = 0. Treating t as a parameter, four solutions for h(y, z, t) are 

h{y,z,t) = ±yd y - 1 g t (y,z,t), (94) 



Hy> z >t) = \d~ l (ygt) = l{ydy 1 g t {y,z,t)-d y : 2 gt{y,z,t)), 
h(y,z,t) = \zd~ l g t (y,z,t), 

h(y,z,t) = \d~ l (zg t ) = ±(zd z 1 g t (y,z,t)-d z 2 g t (y,z,t)). 
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(95) 
(96) 
(97) 



This shows how h can be written in terms of g. For every conservation law corresponding 
to the densities in Table El h could be computed using one of the equations in (I94l)-(l9~7l). 



Conserv ation laws for the K Z equati on have been reported in the literature bv ISharomet 



(119891 ) and ISanders and Wangi (jl997al ). However, substitution of th eir results into (El) 
reveale d inaccuracies. After bringing the mistake to their attention, ISanders and Wang 
(ll997bl ) have since corrected one of their conservation laws to match our result. 



6.3. The Camassa-Holm Equation in 2D 
The (2+l)-dimensional CH equation, 



(u t + ku x - u t2x + 3uu x - 2u x u 2x - uu Zx ) x + u 2y 



0, 



(98) 



for it(x) = u(x, y, t) models water wav es flJohnsonl.l2002l). I t is an extension of t he completely 



integrable 1 D CH equation derived bv lCamassa and Holml (I1993I ). A study bv lGordoa et al. 
(120041 ) concluded that (198]) is not completely integrable. 
Obviously, ([HE]) is a conservation law itself, 



0. 



V t (u x - u 3x ) + V x (ku x + 3uu x - 2u x u 2x - uu 3x ) + T> y (u y ] 
It can be written as a system of evolution equations in y. Indeed, 

-{au t + ku x - u t2x + 3/3uu x - 2u x u 2x - uu 3x ) x 



u. 



V, Vy 



(99) 



(100) 



Note that we introduced auxiliary parameters a and /3 as coefficients of the Ut and uu x 
terms, respectively. The reason for doing so is that the CH equation ( I98p does not have a 
scaling symmetry unless we add scales on the parameters a, j3 and n. Our code guided us in 
finding the following conservation laws with functional coefficients, 

'Dt(fu)+V x (±f(lf3u 2 + ku- \u 2 x -uu 2x -u tx ) + (\fy 2 - \fx)(au t + ku x 
+ 3(3uu x - 2u x u 2x - uu 3x - u t2x )j +V y ((±f'y 2 - ^fx)u y - f'yuj = 

^tyfyuj+V x (^fy{\l3u 2 + nu-\ul-uu 2x -u tx ) + y(\fy 2 -^fx)(au t + 
+ 3/3uu x -2u x u 2x -uu 3x -u t2x )j+Vy(y(lfy 2 -lfx)uy+ (\fx-\f'y 2 )u\ = 0, (102) 

where f(t) is arbitrary and without constraints on the parameters. Thus, if we set a = {3 = 1, 
we have conservation laws for (1981). 



ior 



KU-, 



6.4- The Gardner Equation in 2D 

The (2+l)-dimensional Gardner equation iKonopelchenko and Dubrovskyl ()1984h . 

u t = —\c?u 2 u x + 6/3uu x + u 3x — 3au x d~ 1 u y + 3d~ 1 u 2y , (103) 

for it(x) = u(x, y,t) is a 2 D generalization of 

u t ■■ 



\au 2 u x + <of3uu x + u 3x , 



(104) 
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which is an integrable combination of the KdV and mKdV equations due to Gardner. For 
a = 0, (I103P reduces to the KP equation (fl2j) . For j5 = 0, f ll03|> becomes a modified KP 
equation. Adding a new dependent variable, v = d~ x u y , allows one to remove the integral 
terms and replace fj!03|) by the system 

u y = v x , Vy = \u t - ±u 3x - 2(3uu x + au x v + \a 2 u 2 u x . (105) 

For f)103p . we found two conservation laws with constant coefficients, 

V t (u) + V x (\a 2 u 3 - 3j3u 2 + 3auv - u 2x ) + V y (- (fern 2 + 3vj\ = 0, (106) 

V t (u 2 ) + V x (\a 2 u A - 4{3u 3 + 3au 2 v + 3v 2 + u 2 x - 2uu 2x ) +V V (- u(au 2 + 6v)) = 0(107) 

Using the methodology described for the previous examples in this section, we eventually 
found three conservation laws involving a variable coefficient /(£), 

V t (/u) + V x (/(i« 2 n 3 - 3/3u 2 + 3auv - u 2x ) + f'yv^j 

+ V y (-f(f l au 2 + 3v)-fyu)=0, (108) 

V t (u(fu + &yfj) + V x (/(|a 2 « 4 - 4(3u 3 + 3au 2 v + 3v 2 + u 2 x - 2uu 2x ) 
+ ^yf(koc 2 u 3 - 3(3u 2 + 3auv - u 2x ) + \{2xf + \y 2 f)v) 
+ V y (- fu(au 2 + 6v) - \yf{au 2 + 2v) - l{\y 2 f" + 2xf)u) = 0, (109) 

and 

Vt [{%yf + Pf)u 2 + \{\y 2 f" + xf)u) + V x ((f yf + £/)(§ aV - A(3u 3 

+ 3au 2 v + 3v 2 + u 2 x - 2uu 2x ) + \{\y 2 f" + xf){\a 2 u 3 - 3j3u 2 + 3auv - u 2x ) 
+ \fu x + \y{±y 2 f" + xf)v) +V y (- (%yf + Pf){au 2 + 6v)u 
- \{\y 2 f" + xf){au 2 + 2v) - \y{±y 2 f" + xf")u) = 0. (110) 

Setting f(t) = 1 in (11081 and ffTU9l yields (|TU5|I and ffTUTD . respectively. 

7. Using the Program ConservationLawsMD.m 

Before using ConservationLawsMD.m, all data files provided with the program, as 
well as additional data files created by the user, must be placed into one directory. Next, 
open the Mathematica notebook ConservationLawsMD.nb which contains instructions 
for loading the code. Executing the command ConservationLawsMD[] will open a menu, 
offering the choice of computing conservation laws for a PDE from the menu or from a data 
file prepared by the user. All PDEs listed in the menu have matching data files. An example 
of a data file is shown in Figure [TJ 
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The independent space variables must be x, y, and z. The symbol t must be used for 
time. Dependent variables must be entered as Ui, % = 1, . . . , N, where N is the number of 
dependent variables. In a (l+l)-dimensional case, the dependent variables (in Mathematica 
syntax) are u[l] [x,t], u[2] [x,t], etc. In a (3+l)-dimensional cases, u[l] [x,y,z,t] , 
u[2] [x,y,z,t], etc., where t is always the last argument. 

8. Conclusions 

We have presented an algorithm and a software package, ConservationLawsMD.m, 
to compute conservation laws of nonlinear polynomial PDEs in multiple space dimensions. 

In contrast to the approach taken by researchers working with Maple and Reduce, our 
algorithm uses only tools from calculus, the calculus of variations, linear algebra, and dif- 
ferential geometry. In particular, we do not first compute the determining PDEs for the 
density and the flux components and then attempt to solve these PDEs. Although re- 
stricted to polynomial conservation laws, our constructive method leads to short densities 
(free of divergences and divergence-equivalent terms) and curl-free fluxes. 

The software is easy to use, runs fast, and has been tested for a variety of multi- 
dimensional nonlinear PDEs, demonstrating the versatility of the code. Many of the test 
cases have been added to the menu of the program. In addition, the program allows the 
user to test conservation laws either computed with other methods, obtained from the liter- 
ature, or conjectured after work with the code. The latter is particularly relevant for finding 
conservation laws involving arbitrary functions as shown in Sections \5\ and O 

Currently, ConservationLawsMD.m has two major limitations: (i) the PDE must 
either be an evolution equation or correspond to a system of evolution equations, perhaps 
after an interchange of independent variables or some other transformation; and (ii) the 
program can only generate local polynomial densities and fluxes. However, the testing 
capabilities of ConservationLawsMD.m are more versatile. The code can be used to test 
conservation laws involving smooth functions of the independent variables and the densities 
and fluxes are not restricted to polynomial differential functions. 

Future versions of the code will work with any number of independent variables and 
will cover PDEs that are not of evolution type, e.g., PDEs with mixed derivatives and 
transcendental nonlinearities. 
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(* data file d_kd2d.m *) 
(* Menu item 2-10 *) 

(*** 2D Gardner equation from iKonopelchenko and Dubrovskyl (119841 ) ***) 



eq[l] = D[u[l] [x.y.t] ,y] - D[u[2] [x,y,t] ,x] ; 

eq[2] =D[u[2] [x,y,t] ,y] ~(l/3)*D[u[l] [x.y.t] ,t] + (l/3)*D[u[l] [x,y,t] ,x,3] 
+2*beta*u[l] [x,y,t] *D[u[l] [x,y,t] ,x] -alpha*D[u[l] [x,y,t] ,x]*u[2] [x,y,t] 
-(l/2)*alpha A 2*u[l] [x,y,t] A 2*D[u[l] [x.y.t] ,x] ; 

diffFunctionListlNPUT = {eq[l] ,eq[2] }; 
numDependentVariablesINPUT = 2; 
independentVariableListlNPUT = {x,y}; 

The space variables only; ignore t. 
namelNPUT = " (2+1) -dimensional Gardner equation"; 
notelNPUT = "Any additional information can be put here." ; 

parametersINPUT = {alpha} ; 

All parameters without scaling must be placed in this list. 
weightedParametersINPUT = {beta} ; 

Parameters that should have a scaling factor must be placed in this list. 

userWeightRulesINPUT = {}; 

Optional: the user can choose scales for variables. 
rankRho INPUT = Null; 

Can be changed to a list of values if the user wishes to work with several ranks 

at once. The program runs automatically when such values are given. 
explicitlndependentVariablesInDensitiesINPUT = Null; 

Can be set to 0, 1, 2, ... , specifying the maximum degree (m + n+p) of coefficients 

Q t m x n y p in the density. 
formRhoINPUT = {}; 

The user can give a density to be tested. However, this works only for evolution 

equations in variable t. 

(* end of data file d_kd2d.m *) 

Figure 1: Data file for the 2D Gardner equation in (|103l) . 
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